Quenched dynamics of two-dimensional solitons and vortices in the Gross-Pitaevskii equation 
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We consider a two-dimensional (2D) counterpart of the experiment tiiat led to the creation of quasi- ID bright 
solitons in Bose-Einstein condensates (BECs) [Nature 417, 150-153 (2002)]. We start by identifying the ground 
state of the 2D Gross-Pitaevskii equation for repulsive interactions, with a harmonic-oscillator (HO) trap, and 
with or without an optical lattice (OL). Subsequently, we switch the sign of the interaction to induce interatomic 
attraction and monitor the ensuing dynamics. Regions of the stable self-trapping and catastrophic collapse of 
2D fundamental solitons are identified in the parameter plane of the OL strength and BEC norm. The increase 
of the OL strength expands the persistence domain for the solitons to larger norms. For single-charged solitary 
vortices, in addition to the survival and collapse regimes, an intermediate one is identified, where the vortex 
resists the collapse but loses its structure, transforming into a fundamental soliton. The same setting may also 
be implemented in the context of optical solitons and vortices, using photonic-crystal fibers. 



I. INTRODUCTION 

The last decade has brought about a very substantial amount of research efforts in the physics of atomic Bose-Einstein conden- 
sates (BECs) ll42il43ll . These studies have revealed a wide array of interesting phenomena, not only thanks to the precise control 
over experimental settings and the use of accurate and relatively simple theoretical models, which is a unique peculiarity of this 
area lfl2ll29ll . but also due to direct connections to other areas of physics, including superfluidity, superconductivity, quantum 
and nonlinear optics, and nonlinear wave theory. 

One of the main topics for which these connections have been pursued is the study of the nonlinear dynamics of matter 
waves in BECs. Diverse experimental techniques have been used to produce a broad array of matter-wave excitations. In 
particular, phase engineering has been used to create vortices [i36>, jSTjl and dark solitons H Ull [IH EBl ■ Stirring of the BECs 
has led to the formation of vortices ll24l [32ll and vortex-lattices iB [sl ITtIi . The switch of the scattering length, from positive 
(repulsive) to negative (attractive), via Feshbach resonances, has been used to produce bright matter-wave solitons and soliton 
trains ifTTl [30ll46[ l47ll . These modes have been studied extensively, to the extent that numerous reviews (and even books i29tl ) 
are dedicated to bright solitons 11101, dark solitons ll20ll and vortices 

Hill HI- 

In this work, we aim to revisit a fundamental aspect associated with some of the principal experiments used to produce bright 
solitons, especially those carried out by the Rice group ll46l l47ll . (Note that a similar method was used for the creation of 
solitons in optical fibers in the pioneering works in that field ll22ll35ll .) Precisely, we study the modulational instability (MI) 
ll28ll of the fundamental and vortex soliton after the quench, i.e., after switching the BEC system from repulsive to attractive 
interactions. While this mechanism was extremely efficient in the experimentally elaborated cigar-shaped setting, to the best of 
our knowledge it has not yet been systematically explored in quasi-two-dimensional (2D) pancake-shaped BECs. This is the 
subject of the present work, with emphasis on the formation of solitons and solitary vortices. More specifically, we examine 
the results of the sign switching of the nonlinearity from repulsive to attractive in the 2D BEC, trapped via a combination of a 
harmonic-oscillator (HO), i.e., magnetic, and periodic optical-lattice (OL) potentials IrillOl- Our initial condition, prior to the 
quench, represents either the ground state of the system in the re puls ive-interaction regime, or the dynamically stable excited 
state in the form of a vortex with topological charge 5 = 1 ifisl [19112611 . We note in passing that the vortex is stable in the case of 
the HO trap if its pivot is collocated with a local maximum of the OL potential, yet potentially unstable if is placed at a minimum 

iillll. 

The paper is structured as follows. In section II, we introduce the model and some essential features of the numerical analysis. 
In section III, we report conclusions produced by the simulations for the fundamental solitons. A time-dependent variational 
approximation (VA) for this case is elaborated in section IV. In section V, the case of the solitary vortex is considered. A summary 
of the results and a discussion of possibilities for subsequent work are presented in section VI. 

II. THE MODEL AND COMPUTATIONAL APPROACH 

For sufficiently cold and dilute atomic gases, where the mean-field approximation is well-established, the BEC dynamics 
can be described by the mean-field order parameter (wave function), ^'(r, t). Assuming a strongly anisotropic trap, with the 
transverse (x, y) and longitudinal (z) trapping frequencies chosen so that uj^ — ujy = LiJ± <C oj^, the trapped BEC acquires a 
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nearly planar, ("pancake") shape i29ll42ll43[l . This, in turn, permits one to factorize the wave function, 5* = $(z)i/'(a;, y), where 
is the ground state of the respective HO. Next, averaging the underlying 3D Gross-Pitaevskii equation (GPE) equation 
along the longitudinal direction, z, leads to the following reduced (2D) GPE for the transverse component of the wave function 
(see also Refs. iHHIH): 

ihdti' = ' — Ai' + g2B\i\^i^ + V^^t{x,y)i. (1) 
2m 

Here, A = 9^ + 9^ is the 2D Laplacian, m is the atomic mass, and 1720 = .93D / {V^t^o-z) is an effective 2D coupling constant, 
where g^o = Airfi^as/m (Ug is the scattering length), and the longitudinal trapping length is = y/h/muj^. The potential 
14xt {x, y) in GPE ([T]) is a combination of the HO component and a 2D square-shaped OL: 

1 

2' 



VcKt{x,y) = -muj\r'^ - e[cofp-{kx) +cos^{ky)] 



= V^o{r) + VoL{x,y). (2) 

Here, r is the radial variable, and the OL is characterized by its depth Vo and periodicity aql = 7r/fc. The wavenumber 
k = (27r/A) sin {6/2) (i.e., aoL = A/ (2 sin {6/2))), in turn, is controlled by the wavelength of the interfering beams that create 
the lattice and angle 6 between them. 

Measuring length in units ofaoi^/ir = 1/fc, time in units ofwoL = h/Eoh, and energy in units ofi?oL = 2£'rcc = ^^/"^^ol 
(where E^cc is the lattice recoil energy), GPE ([T]) is cast into the following dimensionless form: 

iut + - {uxx + Uyy) + .g|upM - V^u = 0, (3) 



with potential 



V^^l-n'^r'^ - e [cos 2x + cos 2y] . (4) 



In this normahzed GPE, the wave function is rescaled as t/j ^ \/\g2u\/ Eohi/j exp [i{Vo/Eo'L)t], and the sign parameter is 
g = sign((7oL) = ±1, with g = — 1 and g = +1 corresponding respectively, to repulsive and attractive interatomic interactions. 
Further, the lattice depth in Eq. (|4|i is measured in units of 4ii^ioc, while the normalized HO strength is O = aoL/*^! = '^-L /'^ol, 
where a±_ = h/mujj_ is the transverse trapping length. 

We are interested in the dynamics of fundamental solitons and solitary vortices when the nonlinearity is switched from defo- 
cusing ig < 0) to focusing (g > 0). More precisely, we first fix g = — 1 and solve the imaginary-time GPE version jTol] to find 
the respective steady state for a given configuration, and then solve the GPE with g = 1, using that steady state as the initial 
condition. All the parameters stay the same except that g changes sign. In all the presented examples, ft = 0.1 is assumed in 
Eq. yet our findings should be relevant to a wide range of i7's within the pancake setting. 

Another way (a faster one) of finding such steady state that we employed is to plug ansatz u = e~'''*v(a;, y) into Eq. Q 
to derive a nonlinear eigenvalue problem, which is then solved with the Newton's method. Solutions produced by these two 
approaches are found to be virtually identical, with the maximum pointwise difference being between 10^^ and 10^^" for most 
configurations. 

Unless specified otherwise, the Newton's solution is used as the initial condition in all the simulations. Further, the fourth- 
order split-step Fourier method ll6l [39tl is used to solve the GPE in time. The corresponding domain size is [— Stt, 8tt] x [— Stt, Stt], 
with 256 Fourier modes in each direction and time step At = 0.001. The simulations were run up to T = 2000, which is large 
enough to observe the stability of the final states (if they are stable). Figure[T]shows some examples of the initial conditions, i.e., 
the steady states obtained by the Newton's method. 

Control parameters will be the OL strength e and the norm of the initial condition (i.e., the normalized number of atoms) 

^ = J J Wo{x,y)fdxdy. (5) 

The following main features of the solution will be computed: amplitude (1^1)^^^^^^,, the final norm, which may be slightly different 
from the initial one due to radiation losses, and the total angular momentum, 

M ^ {-i) f f ( ^u* ~ ^u) dxdy, (6) 



^06 d6 

where 6 is the angular coordinate. Note that M is conserved in the isotropic system that does not include the OL, and is not 
conserved with the presence of the OL. In particular, for isotropic solutions with integer vorticity S, 

u = e-'''*+*'^*C/(r), (7) 

the relation between the angular momentum and norm is M ~ 2SN. We consider both the fundamental soliton with S = and 
vortex soliton with topological charge S ~ 1, with or without the OL. 
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FIG. 1: (Color online) A typical example of the output produced by the Newton's method. Such states will be used as initial conditions in 
sections III and V. Left panels: 5 = (fundamental solitons); right panels: 5=1 (vortices). Top and bottom panels correspond to e = 0.5 
and e — —0.5, respectively. 



in. GROUND-STATE QUENCH DYNAMICS 

We start with the fundamental state with no topological charge. We perform the nonlinearity-sign switch (quench) for a wide 
range of OL strengths and initial norms. The former is used as a representative parameter associated with the potential, while 
the latter is employed for scanning through the set of initial data. The resulting two-dimensional map of the stability region 
of the ground state is plotted in Fig. [2] There are two regions with one denoted by open squares, wherein the ground state 
persists in the attractive regime in the form of a stable bright 2D soliton (see, e.g., Refs. ||^|23|] for detailed discussions of the 
stability of such states), and the other denoted by filled circles, which corresponds to the catastrophic wave collapse, a well- 
known phenomenon for equations of the nonlinear-Schrodinger type iH |4^. Clearly, the OL plays a critical role towards the 
stabilization the fundamental soliton, since the critical norm increases as the OL gets stronger. When the OL is absent (e = 0), 
the soliton exhibits breathing dynamics when N < 5.81, and will collapse at > 5.91. Note that the collapse threshold 
corresponding to the Townes soliton is A^cr = 5.85 lH]. 

In the presence of the OL, the dynamics is more complex, and in particular, the angular momentum will be generated when 
e < 0. In such cases, the soliton's center initially coincides with a local maximum of the OL potential, hence it will slide down 
from this position. In fact, in the beginning of the simulation, the soliton moves back and forth between its initial position and 
the centers of other OL cells that are located along a straight line. (It is relevant to note that a 2D soliton can travel more than 
one cell, especially when the lattice is not very strong 114511 .) Then the soliton starts to deviate from this straight line at certain 
time, and the generation of the angular momentum commences. The subsequent motion does not follow any simple pattern. The 
trajectory of the soliton's center for one such case is shown in Fig. [3] 
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FIG. 2: (Color online) The stability diagram for the fundamental soliton resulting from the quench of the ground state of the repulsive BEC. The 
squares denote stable configurations that support breathing dynamics, while the dots denote unstable configurations that lead to the collapse. 
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FIG. 3: (Color online) The trajectory of the soliton's center. The OL strength is e — —0.5, while the soliton's parameters are /i — 0.06, A'' = 
6.0294. The left and right panels are shown from different angles to better represent the rather complex motion. 



IV. THE VARIATIONAL APPROXIMATION 



To develop an understanding of the breathing regime exhibited by the simulations, it is reasonable to apply the VA f33^. The 
starting point is the nonstationary GPE (|3]l, but with a time-dependent nonlinear! ty coefficient g{t). This equation can be derived 
from the Lagrangian, L = dx dyC{u), with density 
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+ ^1^*1 + |- 2 + ^ [cos(2x) + cos(2y)] } \u\ 

The variational ansatz for the fundamental state is based on the isotropic Gaussian, 

1 



M(r, t) — A{t) cxp 



,.2 



(8) 



(9) 



where A, W, b and cj) are, respectively, the amplitude, width, chirp and overall phase, which are assumed to be real functions 
of time. Following the standard procedure, we insert the ansatz into density ([8]) and calculate the corresponding effective 
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Lagrangian, 



The result of the calculation is 



oo 



Leff = 27r / Crdr. (10) 







dS 1 ^db N 1 o , 
(it 2 dt 2W^ 2 

^ n^NW^ + 2eNe-^ , (11) 



47r ^2 2 

where the overdot stands for the time derivative, and the norm of the ansatz is, cf. Eq. Q 

r' + oo /' + 00 

1 2 /i2tt/2 



iV= / / dy|w(x,2/) 1^ = 7rA^M/\ (12) 

J —oo J —oo 

The variational equation 6Lctf/S(f) = reproduces the conservation of the norm, dN/dt = 0, hence N may be treated as a 
constant. Then, equation SLcff/Sb = yields an expression for the chirp: 

and the final equation, SLcff/S (W^) = 0, leads to a closed-form evolution equation for the width, in which b is eliminated by 
means of Eq. (fTSl i: 

'^ = '-^^^-n^W-,eWe-^\ (14) 
dt^ 2nW^ 

For constant g > 0, Eq. ( fT4b yields the well-known variational approximation for the critical norm, A^i™"^^ = 2TT/g ifTill . 

Aiming to compare results of the VA to those of the full simulations, we have to note that, as the initial profiles in the direct 
simulations were taken as per the steady state in the case with repulsive interactions, they are much wider than the OL cell. For 
this reason, the model including the OL cannot be adequately tackled by means of the (radially symmetric) ansatz (|9]), which 
does not include the density modulation induced by the OL. Therefore, the comparison is only carried out for the fundamental 
state in the absence of the OL. 

A number of examples of the breathing regime of the resulting fundamental soli tons are displayed in Fig. ID These are obtained 
by solving Eq. ( fT4] i. for which two initial conditions are needed. The first initial condition W{t = 0) is simply the width of the 
initial soliton, while the second initial condition (t = 0) is approximated by a first order finite difference based on the widths 
of the solution in direction simulation at t = and t = 0.001. 

In Fig.m one observes that, for sufficiently small N, the resulting amplitude of the solution closely follows the VA prediction. 
However, as increases, there arises a beating effect in the full GPE dynamics that is presumably not accounted for by the 
simplified VA dynamics of Eq. ( fT4] i. Nevertheless, the VA still captures the principal features of the oscillatory dynamics of the 
width of the fundamental soliton. 



V. THE QUENCH-INDUCED DYNAMICS OF VORTICES 

With the introduction of initial vorticity, the most interesting observation is the rather delicate character of the stability of 
trapped solitary vortices under the self-attractive nonlinearity (see Ref. Issll and references therein). For the vortical initial 
inputs, we have performed an extensive numerical analysis similar to that reported in section III for the simpler case of the 
fundamental state. The respective two-parameter (the OL strength and initial norm) stability region is shown in Fig. [5] Here, 
we categorize both coherent and less coherent (see the examples below) states generated by the vortical inputs as stable when 
the collapse does not occur by T = 2000 (our reporting horizon). When e = (i.e., the OL is absent), the critical norm for the 
quench process is found to be 

N^J=^"> « 11.81, (15) 

(S—l) 

which is essentially larger than the known value, Ni^ax — 7.79, i.e., the boundary of the existence of (numerically) exact 
stable trapped vortices with topological charge 5* = 1 |38]. Thus, the effective stability range of dynamical (breathing) vortices 
may be essentially broader than that of their static counterparts. For the simulation time exceeding T ~ 2000, the critical norm 
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FIG. 4: (Color online) Comparison of the amplitude of the breathing fundamental soliton, generated by the quench, between the VA and direct 
simulations. The left panel refers to TV = 2.14, and the right one to a vicinity of the critical point, for TV = 5.81. 
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FIG. 5: (Color online) The stability region for the vortical initial condition with S = 1. Here, three regimes are identified. Circles represent 
the case when the resulting state is a vortex with charge S = 1. Triangles correspond to an intermediate regime of non-collapsing solutions, 
which, however, do not keep the vorticity. Finally, dots represent collapsing solutions. 



may be found to vary slightly, as the stable solutions found at TV very close to TVd are still in the process of splittings 
and recombinations at T = 2000. These "hesitating" solutions do not constitute a regular vortex soliton, but rather exhibit an 
additional breathing process. 

The dynamics becomes much more complicated when the optical lattice is present. For these cases, the main conclusions are 
as follows. 

(i) The stability region becomes smaller for a weak OL (e = ±0.1). Yet, the stable solutions still behave like a coherent vortex 
ring, provided that norm TV is not too large. With an increase of TV, the solutions mimic the splitting-recombination scenario 
observed in the absence of the OL (see above). Two typical solutions are displayed in Fig.|6] 

(ii) In the case of the OL with a moderate strength (e = ±0.5), there are two regions in which the solution stays stable for up 
to r = 2000. In particular, it is stable when its initial norm lies in two intervals, 

< TV < 8.04, (16) 
17.59 < TV < 25.17. (17) 

In the stability region ( fT6b . there are two kinds of solutions observed. When the norm is small enough, all the solutions are similar 
to the one shown in the top panel of Fig.|2] Basically, it is built of eight peaks forming two groups. The first group includes the 
peaks at the 3, 6, 9 and 12 o'clock positions, and this set is always present. The second group consists of the remaining peaks 
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FIG. 6: (Color online) The dynamics of a typical stable vortex soliton with S = 1 in the presence of a weak OL. Top: e = 0.1, A'' = 4.6673; 
bottom: e = -0.1, iV = 3.4017. 



set along the diagonals, which breathe as a function of time (thus becoming more or less noticeable). The overall topological 
charge is still preserved js^l . The second kind of dynamics eventually leads to the loss of the vorticity, and transformation of 
the vortex into a fundamental soliton, as shown in the bottom panel of Fig. [7] Its location may vary, depending on the initial 
norm of the solution. The largest observed norm contained in such a stable peak is iV = 5.43. When N is increased to be out 
of the stability region (fTSl i. the solution eventually collapses. In the stability region ([TT] ). the solutions for all the considered 
configurations feature four peaks. The largest total norm is = 21.53, with 5.38 in each individual peak. A solution of this 
type is displayed at Fig.|8] However, only some of the stable solutions are true vortices (marked by circles in Fig.|5l), bearing the 
topological charge of 5 = 1. For other solutions, the loss of the global coherence among the peaks occurs in the course of the 
evolution, at the same moment of time when the symmetry-breaking occurs, i.e., the peaks start to have different height. Two 
examples are presented in Fig. |9] with one preserving the vorticity (similar to the vortices on discrete lattice, although there is 
complete lack of the isotropy in the system |25, 34]), and one becoming incoherent and thus shedding the vorticity off. 

(iii) For large strength of the OL (e.g., e = ±0.9), the first stability region, corresponding to Eq. (fT6] l, expands, similar to 
what is the case for the fundamental states in Section III. On the other hand, the second stability region, which corresponds to 
Eq. (1% . practically disappears. The final stable solutions (at T ~ 2000) exhibit breathing behavior. Yet aside from their more 
pronounced peaks, they do not exhibit any salient structural differences from their counterparts considered above at intermediate 
values of e. 



VI. CONCLUSIONS AND FUTURE CHALLENGES 



In this work, we have examined the quenched dynamics of both the fundamental states and vortices with topological charge 
5 = 1 in BEC. The quench consists of the sudden reversal of the nonlinearity sign from repulsive to attractive. The resulting 
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FIG. 7; (Color online) A vortex soliton with topological charge S — 1, for the OL strength e = 0.5. Top: A*' — 1.3383, the final solution still 
being a vortex. Bottom: = 3.9573, the solution evolving towards a fundamental soliton. 




FIG. 8: (Color online) A vortex solution from the stability region l ll7t . Here the lattice strength is £ = 0.5. 
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FIG. 9: (Color online) Examples of the loss of symmetry for the quasi-discrete vortices from the stability region l ll7t . The OL strength is 
£ = 0.5. Left panels: the largest difference among the four peak amplitudes. Right panels: phase shifts (in units of tt) between the peaks. Top: 
iV = 17.92; bottom: iV = 25.17. 



states were investigated by systematic simulations, addressing both the impact of parameters, such as the strength of the OL 
(optical lattice), and the effect of initial conditions, by considering a range of values of the initial norm, N (which is proportional 
to the number of atoms in the BEC). A principal result is that the OL expands the range of initial norms which do not lead to the 
collapse of both the fundamental and vortex states. In fact, this expansion occurs well over the interval of norms for which static 
vortices were previously identified as stable states via the linear stability analysis. For the fundamental states, the application 
of the VA (variational approximation) is more efficient for the lower norm of the initial state. Additional beating effects, not 
captured by the VA, were found close to the collapse threshold. On the other hand, a particular finding, in the case of the vortex, 
is that, in addition to the regimes where a vortex survives in the breathing form or collapses, there is an intermediate regime 
where the vortex (in the presence of the OL) loses its topological character, yet remains immune to the collapse. Furthermore, a 
second stability region was identified for the vortex, in the presence of the OL with an intermediate strength, where the collapse 
was avoided due to the formation of a robust quasi-discrete vortex. 

The same system may also be implemented in nonlinear optics, where the combination of the HO trap and OL potential 
corresponds to photonic-crystal fibers, while the switch of the nonlinearity corresponds to a junction of two waveguides made 
of self-defocusing and self-focusing materials ll49ll . 

We believe that these results provide a potential for further studies on this theme, both at the theoretical and at the experimental 
level. At the theoretical level, it may be useful to develop 3D generalizations of the analysis developed in this paper, either for the 
full case of spherically symmetric BECs or for strongly anisotropic traps. The effect of anisotropy in 2D would be interesting 
to examine too, as it departs from the configuration of the isotropic cylindrical trap. On the other hand, at the experimental 
level, recent advances in producing lliHlsoll and monitoring ll2lll37ll the dynamics of vortices and vortex clusters, in conjunction 



with the well-established control of the BEC dynamics by means the Feshbach resonance |44ll , render particularly appealing 
and accessible examination of such quenches (or the corresponding adiabatic transitions) between the repulsive and attractive 



10 



regimes. 
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